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Full charge self-consistence (CSC) over the electron density has been implemented into the local 
density approximation plus dynamical mean-field theory (LDA+DMFT) scheme based on a full- 
potential linear muffin-tin orbital method (FP-LMTO). Computational details on the construction of 
the electron density from the density matrix are provided. The method is tested on the prototypical 
charge-transfer insulator NiO using a simple static Hartree-Fock approximation as impurity solver. 
The spectral and ground state properties of bcc Fe are then addressed, by means of the spin- 
polarized T- matrix fluctuation exchange solver (SPTF). Finally the permanent magnet SmCos is 
studied using multiple impurity solvers, SPTF and Hubbard I, as the strength of the local Coulomb 
interaction on the Sm and Co sites are drastically different. The developed CSC-DMFT method 
is shown to in general improve on materials properties like magnetic moments, electronic structure 
and the materials density. 



I. INTRODUCTION 

Understanding the interplay between strong electron- 
electron interactions and kinematic effects in crystals 
pose a major experimental and theoretical challenge. Al- 
though the problem is very demanding, the potential re- 
wards are great. Many interesting materials properties, 
like high temperature superconductivity[l[ and coexis- 
tence of superconductivity and magnetism, are found in 
materials with a correlated electronic structure, e.g. the 
newly discovered Fe-pnictidesQ. Strong correlations are 
also associated with mixed valence and heavy Fermion 
materials. 

On the theoretical side it now stands clear that a 
straightforward application of conventional electronic 
structure theory 0, as given by the local density ap- 
proximation (LDA) or generalized gradient approxima- 
tion (GGA), often fails in describing correlated electron 
materials. In addition it is evident that an approach 
based on model Hamiltonians also has drawbacks in that 
many of the parameters of these models are hard to es- 
timate and that the band-formation often is described 
inaccurately. Hence it becomes important to combine 
the best features of the two approaches, something which 
can be achieved by combining density functional theory 
in the local density approximation and dynamical mean 
field theory (LDA+DMFT) @, @|. To date there have 
been several reports that integrate DMFT with exist- 
ing codes for electronic structure calculation. Starting 
with the pioneering papers of Ref@4I3j various imple- 
mentations have been produced, differing in the parent 
electronic structure method and the way the local or- 
bitals are constructed. At the beginning most of such 



implementations were limited to work with a fixed LDA 
electron density. This was not only due to the tech- 
nical difficulties in constructing a reliable electron den- 
sity from the DMFT Green's functions, but also to the 
fact that fully charge self-consistent (CSC) simulations 
would have needed a prohibitive amount of computa- 
tional resources. However, in the last years these lim- 
itations have been overtaken by following two different 
approaches. A straight forward but approximate com- 
putational scheme has been proposed by associating the 
new electron density to an LDA+U problem [ll| con- 
structed with the DMFT density matrix [II, EDI . Several 
fully self-consistent LDA+DMFT suites have also been 
produced [14H18I ] , although it remains to be seen if their 
applicability is limited by the high computational cost. 

In the present paper we report the extension of an 
LDA+DMFT implementation[ig|-l2l[ based on the full- 
potential linear muffin-tin orbitals (LMTO) methodp3- 
|25] in which we include full self-consistency in the elec- 
tronic density. Due to the full potential character such 
an implementation is suitable for simulations of any sys- 
tem, without restriction on the geometry. In addition 
the usage of a small number of basis functions associated 
to the LMTO method makes it accessible to calculate 
the fully hybridizing electronic structure of very large 
systems without any ad-hoc assumptions. The technical 
details of the current development, together with a brief 
review of the DMFT equations, are presented in Section 
II, and in Appendix A. 

In section III, we apply the developed LDA+DMFT 
scheme to several materials with various degrees of 
electron correlation, to illustrate the usefulness of the 
method. Our first test-case is the prototypical insulator 



NiO in its sodium-chloride phase. This is an archetypi- 
cal compound which is often used to illustrate that elec- 
tronic structures based on LDA or GGA alone fail[26l|. 
Our second test case is represented by the itinerant fer- 
romagnet bcc Fc. Here the effects of electron correlations 
are much less pronounced, but several well known failures 
with conventional LDA/GGA theory are known, e.g. the 
too broad band widths of the 3d states 27| and the pres- 
ence of a weak intermediate-energy satellite 11 9|. which 
is observed in photoemission experiments [28]. Finally 
we study the electronic structure of the rare-earth based 
permanent magnet SmCo5 , which is a particularly useful 
case, since the correlated electrons on the Sm and Co sites 
experience different strenght of the local Coulomb inter- 
action, and therefore require different impurity solvers. 

In the last section of this paper, section IV, we present 
our conclusions and perspectives of further development 
and future applications. 



II. THE CHARGE SELF-CONSISTENT 
LDA+DMFT SCHEME 

The starting point of the LDA+DMFT scheme is the 
assumption that the single-particle LDA Hamiltonian 
jjLDA foes not correctly describe the exchange and cor- 
relation effects associated to a certain set of local orbitals 
{\R, £)}, where R denotes the lattice site and £ the quan- 
tum numbers labeling such states. A local two-particle 
operator Ur, describing an effective Coulomb interac- 
tion between the Kohn-Sham quasi-particles, is therefore 
added to H LDA , giving the so-called LDA+U Hamilto- 
nian [nisiiii: 
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FIG. 1. The charge self-consistent LDA+DMFT scheme. The 
LDA+DMFT cycle is the union of the LDA cycle (lower sec- 
tion) and the DMFT cycle (upper section). A charge non- 
self-consistent (one-shot) LDA+DMFT scheme lacks the con- 
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electronic bath. This leads to the LDA+DMFT scheme 
depicted in Figure Q] It consists of two interconnected 
self-consistency cycles, the LDA cycle and the DMFT 
cycle, which are presented in more detail below. 



From this point on the site subscript R denotes that the 
(eventually k-dependent) operator is projected onto the 
local orbitals 



a r = y, \R, (R, CI E A *\ R > 0( R > ?\ 



(2) 



The double-counting correction H^ c in Equation (fT]) is 
introduced to remove the improper LDA description of 
the local Coulomb interaction. It is clear that there are 
many issues related to the arbitrariness of the choice of 
the correlated orbitals, and to the difficulty of determin- 
ing a reliable double-counting correction. However, a 
discussion of such topics would be out of the scope of 
this paper, so we refer the interested reader to some re- 
cent reviews [Ell] for more details regarding the LDA+U 
Hamiltonian and the DMFT. 

The main idea of the DMFT is that the Hubbard 
model described by the LDA+U Hamiltonian can be 
mapped locally to an effective Anderson impurity model 
(AIM)Q. In this way the lattice problem is reduced to 
the much simpler problem of an atom embedded in a 



A. The LDA cycle 

The starting point of the LDA cycle is the calculation 
of the effective LDA potential V LDA from some electron 
density p(r), using the LDA functional. V LDA is then 
used to construct a k-point dependent LDA Hamiltonian 
jjLDA rjv^g p ure LDA cycle can then be closed by con- 
structing a density matrix 



Pk 



LDA = i + exp {-l3{Ht DA - pi) 



(3) 



from which a new electron density p(r) can be obtained. 
In Equation ((3]) the chemical potential p is set to ensure 
that the system contains the correct number of electrons 
at temperature T = 1//3. This pure LDA cycle is shown 
in the lower part of FigfTJ After the LDA cycle has pro- 
duced a reasonable initial electronic structure, one may 
exit the LDA cycle and continue into the DMFT part of 
the full LDA+DMFT scheme. Once the calculation has 
entered the DMFT cycle, the continued use of the LDA 
cycle will be limited to updating H^ DA from p DMFT . 
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Details of how to construct the electron density from a 
representation of the density matrix in the LMTO basis 
are presented in Appendix [A] 



problem is to decompose the Green's function into a nu- 
merical part and an analytical part, taking advantage of 
the asymptotic expansion of the self-energy for large w n 



B. The DMFT cycle 

The DMFT cycle starts with the creation of the k- 
dependent Green's function Gk(iu> n ) from Ht DA and a 
set of local self-energies {Yi R {iw n )} 



G k {iuj n ) 



iu} n l 



(4) 

where Gk(iuj n ) is resolved over the Matsubara frequen- 
cies ioj n = inT(2n+l) from the finite temperature many- 
body formalism. At the first iteration the k-independent 
self-energies are set to zero, but in the following ones 
they will be obtained from the solution of the effective 
impurity problem. The chemical potential /i in Equa- 
tion ([4]) needs to be readjusted at every DMFT iteration 
to ensure that the number of electrons in the system is 
conserved. The mapping to the effective impurity model 
is done through the use of local bath Green's functions 
Go R (iu) n ), which are obtained by means of an inverse 
Dyson equation, 



& 0R (iu n ) = G R l {iuj n ) + fl R (iu n ). 



(5) 



The U term in Equation ([T]) and the bath Green's func- 
tion GvR,(i<jo n ) fully describes the effective AIM. This 
model is well studied in literature and there exists many 
impurity solvers. The resulting set of impurity self- 
energies 'E A / M (iuj n ) are then used to update the local 
self-energies after the double counting contribution has 
been removed 
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The DMFT cycle can now be continued using the up- 
dated self-energies, as shown in the upper part of Fig. [T] 
Alternatively, a new density matrix p^ MFT can be calcu- 
lated from the k-dependent Green's function, to update 



H£ DA in the LDA cycle. 



C. Updating the electron density 

The density matrix p^ MFT can be obtained from the 
Green's function as an infinite sum over all the Matsub- 



ara frequencies [31 



N 

p k = lim lim TV [G k (iw„)e !IJ ™" + G[(iw„)e- !1J " 

n— 

The order of the limits can not be interchanged as the 
partial sums are only point wise convergent and not uni- 
formly convergent with respect to r). A way around this 
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the density matrix can be split into two parts as Pk = 
Pk" + Pk um i where 
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The analytical part p£ n has a simple form but contains 
the logarithmic divergence of pk, while p™ um converges 
uniformly (2lj . The uniform convergence allows the or- 
der of the limits in Equation (|13j) to be interchanged. 
With a minimal loss of accuracy the resulting sum can 
be truncated at some large cut-off Matsubara frequency 
Nmax, giving 



N„ 
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Pk 
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(13) 



The analytical part of the density matrix can still not be 
summed explicitly, but thanks to its frequency indepen- 
dent form it can be evaluated in the same way as p\P A 
in Equation ([3]), which yields 



i + exp(^(i*k^ + E^M-4)^ • 

(14) 

Notice that in the absence of a self-energy G^ um (iu) n ) — 
0, and p£™ is reduced to p^ DA - 



D. The asymptotic limit of E_r(iu;„) 

The problem of evaluating the infinite sum in Equation 
([7]) has now been replaced by finding the value of the 
Hartree-Fock-like self-energy S^(co) from the asymp- 
totic limit of Yl R (iuj n ). The asymptotic behavior of 
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S^(io; n ) is dictated by the Dyson equation in Equa- 
tion (|5|). Unless iV max is set to a very large value, at 
a high cost of computational resources, the slow decay of 
the anti-hermitian part of 'Sn(ioj n ) makes the zeroth or- 
der approximation £#(00) ~ Sfl(iwjv max ) far too crude. 
However, the situation can be improved thanks to the 
hermiticity of Er(oo). The contribution from the slowly 
decaying tails can be canceled out by taking 



(15) 

The current implementation goes one step further and 
extrapolates Sj|(iw n ) to infinity from a least squares fit 
of a few data points around iV max . Since Xj|(iu;„) is 
hermitian it can be expanded in orders of w~ 2 , so the 
first order fit has an error proportional to w^- . 



III. RESULTS 

The following calculations were performed with an 
LMTO basis set containing a triple basis for s and p states 
and a double basis for d and f states, where the multiplic- 
ity is referred to the number of envelope functions (tails) . 
The basis of the valence electrons was constructed con- 
sidering 4s, 4p and 3d states for the transition metals 
atoms, 2s, 2p and 3d states for the O atoms, and 6s, 6p, 
5d, and 4f states for the Sm atoms. The correlated basis 
functions were chosen to be the MT orbitals described 
in Ref. [U The k-points were distributed according to 
the conventional Monkhorst-Pack grid, and the Brillouin 
zone integration was carried out using Fermi smearing 
for the temperature T = 400 K. Finally the U-matrix in 
Equation ([1]) was constructed from the Slater parame- 
ters F°, F 2 , F 4 , and F 6 Q- For Fe and Co an additional 
parametrization has been adopted by using fixed atomic 
ratios for F 2 and F 4 , which allows to refer to the bare 
Coulomb repulsion U — F° and to the Hund's exchange 



A. Nickel oxide 

The first test of our CSC LDA+DMFT implemen- 
tation is the study of the antiferromagnetic (AFM) 
NiO within the LDA+U approximation. Although the 
LDA+U scheme does not necessarily require the usage of 
a Green's function formalism, it can still be rewritten in 
terms of the equations presented in Section IIIB1 solving 
the effective impurity model by means of the Hartree- 
Fock approximation. NiO in LDA+U is a good test since 
it is the standard test-case used in most LDA+U imple- 
mentations and much material is available in the liter- 
ature. Moreover the resulting self-energy is static, and 
therefore calculations can be made using only a single 
Matsubara frequency. This implies that the numerical 
contribution p" um is zero by construction. Finally it is 



important to stress that a proper description of the elec- 
tronic properties of NiO would require a more sophisti- 
cated solver, especially if a detailed comparison with ex- 
perimental data, e.g . from photoemission spectroscopy, 
is under focus (34l - l36j |. 

The crystal structure of NiO is NaCl but due to the 
antiferromagnetic structure, the calculations involved 4 
atoms per unit cell, 2 O and 2 Ni. The lattice constant 
used was 7.89 a.u. and the Brillouin zone was sampled 
through a mesh of 9 x 9 x 9 k-points, giving 365 points 
in the irreducible wedge. Spin-orbit coupling was consid- 
ered with spin quantized along the easy axis (111). The 
values F = 8.00 eV, F 2 = 8.19 eV and F 4 = 5.11 cV 
were taken from Ref. H3- Finally the double-counting 
correction was chosen to be the fully localized limit [HJ , 
as is usual for insulators. 




FIG. 2. (Color online) LDA and LDA+U spin-integrated 3d 
density of states of NiO in the present RSPt implementation 
and in the ELK code[3S|. The computational details are de- 
scribed in the main text. Notice that the curves have been 
shifted in order to allign the first occupied peak, due to the 
arbitrariness of the chemical potential in insulators. 



The electronic structure is shown in Figure [2] for both 
the LDA and the LDA+U approximations. The figure 
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TABLE I. Ni 3d contribution to the spin moment, m s , orbital moment, m , and total moment, m to t, for NiO in RSPt and 
ELK. These values are compared with the experimental data from Ref. [3^ and [4(J. All the magnetic moments are given in 
units of /is per atom. 



Method 


RSPt 

m s m D m to t 


ELK 

m s m Q m to t 


LDA 

LDA+U (CSC) 
Exp. [39, 40] 


1.21 0.14 1.35 
1.73 0.24 1.97 
1.7^2.1 0.28^0.37 2.0-^2.4 


1.14 0.14 1.28 
1.70 0.26 1.96 
1.7^2.1 0.28-^0.37 2.0-^2.4 



contains a comparison between results of the present 
implementation and the ELK code [HI, which is based 
on the full-potential linear augmented plane-wave (FP- 
LAPW) method. As is clear from the figure, the elec- 
tronic structures of both methods are practically identi- 
cal when LDA is considered. When LDA+U corrections 
are introduced the main features of the spectral proper- 
ties are also in good agreement, but some smaller differ- 
ences may be found for the finer details. These discrep- 
ancies can be linked to a different definition of correlated 
orbitals, which leads to a slightly stronger effect of the 
Hubbard U for the RSPt results in comparison with the 
ELK calculations. This can be seen in the peak struc- 
ture at -7.5 eV, which is found to lie a bit lower in the 
RSPt implementation, and the unoccupied peak at 3.5 
eV, which is a bit higher. For both implementations the 
electron-electron interaction, parametrized in the form 
of U, widens the band-gap, and the overall agreement of 
these results with previous LMTO-based studies is excel- 
lent dum. 

Some additional information can be gathered by look- 
ing at magnetic moments, shown in Table fl] We first no- 
tice that the agreement between RSPt results and ELK 
calculations is very good, and both methods give results 
compatible with recent experimental data[39l |40|. 



B. bcc Fe 

We turn next our attention to a system where the 
electron-electron interaction among the 3d states is 
weaker than that of NiO, namely bcc Fe. In general for 
the itinerant ferromagnets a static solver as the Hartree- 
Fock method used in the previous section is not sufficient, 
since the energy dependence of the self-energy plays a rel- 
atively important role. Instead we have for these calcu- 
lations used the more sophisticated SPTF solver [HI, [43[ . 
Since correlation effects are weak, we do not expect big 
changes associated to the CSC cycle in the LDA+DMFT 
scheme. However the self-energy is now a dynamical 
quantity, and therefore the numerical contribution to the 
electron density p" um becomes finite. 

The spectral properties of bcc Fe were calculated at 
the experimental lattice constant [l9| and with spin-orbit 
coupling with the spin axis aligned along the easy axis 
(001). The Brillouin zone was sampled through a mesh 
of 21 x 21 x 21 k-points, leading to 726 points in the 



irreducible wedge. Finally the averaged static part of the 
self-energy was used as double-counting correction [32| . 




FIG. 3. (Color online) Projected 3d density of states for bcc 
Fe in the RSPt code for basic LDA, one-shot LDA+DMFT 
and CSC LDA+DMFT for U=1.7 eV. The Fermi level is at 
zero energy. 

In Figure [3] we show the spin-resolved 3d density of 
states of bcc Fe for plain LDA, one-shot LDA+DMFT 
and CSC LDA+DMFT for U=1.7 cV and J=0.9 eV. Such 
a choice of U and J was made considering the recently 
published LDA+DMFT study of the spin- and angular- 
resolved photoemission spectra of Fe (110) [27j, and the 
successful application of the LDA+DMFT scheme to the 
bcc-to-fcc (a-7) phase transition of Fe [44]]. As observed 
in previous works [ill H^, [4(| the local correlation ef- 
fects are stronger for majority spin states, making the 
3d band narrower, which is consistent with experiments 
27]. Moreover the formation of a high energy satellite 
can be noticed, which has been thoroughly discussed in 
Ref. [Til. It is relevant to focus on the fact that essen- 
tially no differences can be observed between one-shot 
and CSC results for bcc Fe. In order to explore the role 
of the CSC cycle, we have made additional simulations 
for higher values of U (not shown here). Even for U 
as high as 4 eV no significant differences can be seen 
between one-shot and CSC densities of states. This is 
related to two factors. First of all our choice of the dou- 
ble counting cancels much of the effects of the Coulomb 
interaction, similarly to the so-called around-mean field 
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TABLE II. 3d contribution to the occupation n, spin moment 
m s , and orbital moment m Q for bcc Fe. RSPt results for 
plain LDA, one-shot LDA+DMFT and CSC LDA+DMFT 
are presented for the two chosen values of U. All the magnetic 
moments are given in units of fis per atom. 



Method 


one-shot 
n m s m Q 


CSC 
n m s m Q 


LDA 

U=1.7 eV 
U=2.3 eV 
U=3.5 eV 
Exp. [47, 48] 


6.24 2.21 0.05 

6.27 2.17 0.07 

6.28 2.19 0.06 
6.32 2.25 0.05 

- 2.13 0.08 


6.24 2.21 0.05 

6.25 2.14 0.06 
6.25 2.15 0.06 
6.27 2.22 0.04 

- 2.13 0.08 



0.12 



- 0.08 



< 



0.04 



0.00 




100 

Atomic Volume (a.u/) 



140 



FIG. 4. (Color online) Total energy versus atomic volume 
curves of bcc Fe for LDA, one-shot LDA+DMFT and CSC 
LDA+DMFT for two different values of the local Coulomb 
interaction . The total energies are plotted with respect to 
their minima. 



double counting for LDA+U (ll|. Additionally the fact 
that the 3d states only hybridize with very flat sp bands, 
makes the electron density insensitive to small changes 
induced by the self-energy. 

Despite the fact that the spectral properties look simi- 
lar for the one-shot and the CSC calculations, differences 
can be observed when looking at the occupation num- 
bers, spin and orbital moments, as show in Table [TTJ The 
3d occupation in one-shot LDA+DMFT tends to increase 
when compared to the bare LDA value, while the CSC 
calculation results in an occupation close to the LDA re- 
sult. Such a behavior can be attributed to the fact that 
the H^ DA is now fixed, and not allowed to adjust to com- 
pensate the changes induced by the local interactions in 
the electron density. Moreover from Table |TT] we can also 
notice how the LDA+DMFT values improve on the LDA 
results for the spin and orbital moments, when compared 
to experimental data, especially for U=1.7 eV. However, 
the agreement is not perfect, and we identify the problem 
in the perturbative nature of the SPTF solver. Starting 
from the SPTF functional other computational schemes 
can be constructed, which do not suffer of such a problem 
491. 



Although the differences of the occupation numbers 
due to the CSC cycle are small, we may expect more 
important changes to be seen in the total energy correc- 
tions to the bare LDA functional. Following the scheme 
presented in Ref. [H we calculated the LDA+DMFT 
total energy for a wide range of lattice constants. To 
reach a convergence up to the meV, a higher number of 
Matsubara frequencies was needed, i.e. 8192 frequencies, 
leading to an energy cut-off of about 20 Ry. The energy 
versus atomic volume curves for LDA and LDA+DMFT 
are shown in Figure HI Already for such a small value of 
the Coulomb interaction as U—1.7 eV the changes on the 
ground state properties of the material are strong: the 
minimum is shifted to higher atomic volumes and the 
curve is flattened. Contrary to what one might expect, 
also for the total energies the differences between one- 
shot and CSC calculations are very small, and the corre- 
sponding curves in Figure 2] are almost coinciding. For 
an ideal description of the ground-state properties of bcc 
Fe within our theory a smaller value of U is needed: the 
data for LDA+DMFT total energies for U=1.2 eV and 
J=0.8 eV are also shown in Figure HI A better overview 
of these results can be obtained by looking at Table IIII1 
where the equilibrium atomic volume and the bulk modu- 
lus are reported. Such values were obtained by fitting the 
total energy data through the Birch-Murnaghan equa- 
tion of state [13, [H|. Fitting through the Vinet equation 
of statefHJ was also made, leading to very similar results 
(not shown here) . As for our previous study of fee Ni [2l[ , 
it is clear that our calculations tend to over-estimate the 
role of the local exchange and correlation effects in the 
total energy. We identify two possible reasons for such 
a behavior. The first one can be found in the perturba- 
tive nature of the SPTF solver, which is supposed to be 
used for systems where the Hubbard U is smaller than 
the bandwidth W (42j. When increasing the lattice con- 
stant the bandwidth decreases, and the applicability of 
our approximations becomes questionable, at least con- 
cerning such a subtle feature as the total energy. The 
second reason is that the SPTF solver is usually applied 
in a non-conserving form, in Baym-Kadanoff sense [49j , 
and this can introduce an additional error in the total 
energy. Although such a feature had already been ob- 
served for fee Ni, it is especially interesting to see that it 
is present also for bcc Fe, where the bandwidth is bigger 
and the Coulomb repulsions are weaker. 



C. SmCo 5 

S121C05 has the CaCus type crystal structure given 
by space group #191 (P6/mmm) with lattice parame- 
ter a = 5.002 A and c = 3.964 A, with Sm in Wyckoff 
position la, and Co both in 2c and 3g. The k-point mesh 
was set to 8 x 8 x 12, giving 448 points in the irreducible 
wedge of the Brillioun zone. 

In the LDA density of states, as shown in Fig. \5jp, the 
Sm 4f states are found in two narrow spin split peaks at 
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LDA 
LMTO LAPW 


U=1.2 eV 

one-shot 


U=1.7 eV 
one-shot CSC 


Exp. 


Vo (a.u/) 
B (GPa) 


70.49 70.71 
234 233 


78.51 
148 


87.00 87.06 
91 90 


79.51 
168 



TABLE III. Equilibrium atomic volume Vo and bulk modulus 
B from LDA, one-shot LDA+DMFT and CSC LDA+DMFT 
for two different values of U. These values were obtained by fit- 
ting the calculated data through the Birch-Murnaghan equa- 
tion of stat e [50l . l5l| for a temperature T = 400K. Experi- 
mental data[53l] and results of previous full potential LAPW 
studies 54] are shown for comparison. 



and 3.4 eV. The small bandwidth indicates that the f- 
electrons are strongly localized and should be treated as 
atomic-like. The Hubbard-I Approximation [kJ HH HI] 
(HIA) has been shown to give excellent results for 
strongly correlated f-orbitalsjsilj [E3, HH and was there- 
fore chosen as the impurity solver for the Sm 4f orbitals. 
The Co 3d states form a broad peak between -4 to eV in 
the LDA density of states, close to what is found in Co 
metal [l9L [59j . which points to weak correlation effects. 
With this in mind, the SPTF impurity solver in the fully 
renormalized version presented in Ref. l49l was chosen to 
treat the Co 3d states. 

The double counting correction H^ c — fJ-Dci + Hx 
was used for the Sm 4f states, where the double count- 
ing parameter /idc acts as an atomic chemical potential, 
and Hx was constructed to remove the local intra-orbital 
LDA exchange splitting. It is not possible to fully deter- 
mine fiDC from the number of Sm f-electrons alone, un- 
less the system is intermediate valence 20]. In the present 
calculation the value of [luc was set by requiring that 
the Sm ground state should contain 5 electrons, and that 
the lowest binding energy of the / 5 configuration should 
be at 5.7 eV. Hx was constructed under the assump- 
tion that the Sm 4f inter-orbital LDA exchange splitting 
was only slightly perturbed by the polarization of the f- 
states. An estimate of the inter-orbital LDA exchange 
splitting could then be obtained from a separate charge 
self-consistent calculation, but where the Sm 4f electrons 
were constrained to stay paramagnetic. One-shot calcu- 
lations would completely neglect the effect of the SPTF 
self-energy on inter-orbital exchange splitting, and there- 
fore charge self-consistence is necessary for an accurate 
calculation. The double-counting for the Co atoms is 
chosen to be the same as that for bec Fe in the previous 
section. 

The values of the Slater parameters F 2 , F 4 and F 6 
for the Sm 4f orbitals were recalculated at each LDA 
iteration from the radial integration of the bare Coulomb 
interaction over the f-orbitals. To account for some of the 
screening by the non-f-electrons, the values of the Slater 
parameters F 2 and F 4 were reduced to 92% and 97% of 
their calculated values, respectively. A reliable value of 
the Slater parameter F° can not be obtained in this way, 
as it is strongly affected by the screening. Therefore F° 



TABLE IV. Summary of the impurity solvers (AIS) and 
final Slater parameters (F°,F 2 , F 4 and F 6 ) used in the 
LDA+DMFT calculation of SmCo 5 . 



Corr. orbitals 


AIS 


F u (eV) 


F^ (eV) 


F 4 (eV) 


F° (eV) 


Sm 4f 


HIA 


8.00 


11.63 


7.64 


5.65 


Co 3d 


SPTF 


2.50 


7.75 


4.85 





TABLE V. Spin (m s ), orbital (m G ), and total (m tot ) moment 
of SmCos compared with experimental data at T=0 K. All the 
magnetic moments are given in units of \ib per atom, while 
the final total moment is per formula unit. The contribution 
from the interstitial regions amount to 0.4 hb- 





Sm 


Co (2c) 


Co (3g) 


Total 


Method 


m s m G 


m s m D 


m s m D 


m to t 


LDA 


-5.48 1.81 


1.58 0.06 


1.55 0.10 


4.04 


LDA+DMFT 


-3.47 3.26 


1.54 0.22 


1.52 0.18 


8.02 


Exp. [62] 








8.9 



for the Sm 4f state was set to 8 eV, following 
The Slater parameters for the Co 3d orbitals were taken 
from previous LDA+DMFT studies jl9|] . A summary 
of the solvers, and the Coulomb parameters used in the 
LDA+DMFT calculation is given in Table [TV] 

The projected density of states of the Sm 4f and Co 3d 
orbitals from the CSC LDA+DMFT and LDA calcula- 
tions are shown in Figured Both LDA and LDA+DMFT 
correctly position the Co 3d states at the Fermi level, but 
only the latter gets the shape of the peak in excellent 
agreement with the XPS spectrum of Ref. The dif- 
ferences are much larger for the Sm 4f states, where LDA 
yields two narrow spin-split peaks instead of the multi- 
plet structures found in LDA+DMFT. The position of 
the first Sm 4f multiplet peak in the LDA+DMFT cal- 
culation at -5.7 eV coincides with the peak of the XPS 
spectrum, since the double counting parameter [Ibc was 
tuned with respect to this peakposition. The lack of 
structure between -10 to -7 eV in the experimental XPS 
spectrum is puzzling, as one would expect to see signs of 
the remaining multiplet structure like in Sm metal and 
SmCo 2 [6l|. 

Finally, the magnetic moments are shown in Table V 
and the total moment from the theory is 10% lower than 
the experimental data at zero Kelvin. Theoretical val- 
ues obtained by LDA are also shown, but they are much 
smaller than in experiment. The data in Table V are re- 
warding in that for the Sm atom, the Russel-Saunders 
limit seems to be approaching, which is expected for a 
reasonable theoretical treatment of rare-earth based com- 
pounds. 

IV. SUMMARY AND CONCLUSIONS 

We have developped a charge self-consistent 
LDA+DMFT code, and have applied it to the study of 
the electronic properties of a few weakly and strongly 
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-10 -5 5 10 15 

Energy (eV) 



FIG. 5. (Color online) Experimental and calculated spectra 
of SmCo5 , with the Fermi level at zero energy, a) Experimen- 
tal x-ray photoemission spectrum from Ref. [6C|. The small 
peak at -11 eV is an Ar 2p artifact from the sample cleaning 
procedure, b) Calculated projected density of states of Sm 
4f orbitals in LDA (dashed red line) and LDA+DMFT (thick 
black line), c) Calculated projected density of states of Co 3d 
orbitals at the Wyckoff position 3g in LDA (dashed red line) 
and LDA+DMFT (thick black line), d) Calculated projected 
density of states of Co 3d orbitals at the Wyckoff position 2c 
in LDA (dashed red line) and LDA+DMFT (thick black line). 



needed to define a proper double-counting. Excellent 
agreement was found with existing photoemission data 
for low binding energies, but not for higher energies. In 
perspective we plan to apply the present methodology to 
study more complex properties of permanent magnets, 
complex oxides and corelated f-electron systems. Such 
studies are underway. 
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Appendix A: The full potential electron density 

In the following Appendix, we present some techni- 
cal details of the calculation of the electron density from 
a density matrix in the FP-LMTO basis. The FP- 
LMTO basis is technically involved, and rather than go- 
ing through all details, we will simply state the relevant 
expressions in the notation of Ref. |23| with only a brief 
explanation of the symbols. 

In a normal single particle LDA calculation the den- 
sity matrix of Equation is obtained by squaring the 
eigenvectors, weighing each vectors contribution by the 
appropriate Fermi weight, 

(i\j% DA \j) = Pij (k) = Y, w ".'« A *(vMA}(v,bk). (Al) 

V 

Here Ai(v, k) is basis function i component of the eigen- 
vector v, and uyk is the corresponding Fermi weight. 



1. Within the muffin-tin spheres 



correlated materials. Antiferromagnetic NiO was used 
as a first test of our implementation, and comparison 
with other codes for electronic structure lead to a good 
agreeement, taking into account the technical differences 
involved in the calculations. Then we have addressed the 
problem of the role of the charge self-consistence in the 
LDA+DMFT study of the weakly correlated bcc Fe. It 
was found that the update of the LDA electronic density 
leads to minute changes of the spectral and magnetic 
properties. The total energy and the bulk modulus 
are also quite unchanged. However for the strongly 
correlated SmCos the charge self-consistence is strictly 



From a density matrix such as that given in Equation 
(IA1[) . the electron density rk(r) inside the muffin-tin of 
type t, is given as the following sum, 

u t (r)=^D w (f)^U^(e',r) 

h e>t' (A2) 

J2M M (e,£;e',£')\J^(e,r), 

el 

where U is a radial basis function and D^t is a symmetric 
linear combination of spherical harmonics. We define the 
two intermediate density matrices, 
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M M {e,£; e',£' 



24 + 1 

47T 



E 

m,m' ,mh 



Q(£, m; £', m';£ h , m h )a* (m h )M t (e, £, m\ e', m') 



(A3) 



Mt(e,£,m;e',£',m') = ^t-tta 6(e, ej)Q t ,l{e, /ti)Si m; i imf (e; k^; t - T- j; ;k)x 

p y (k)n t ,i/(e', «j)S], m/ .^ mj (e'; K j; x - T,-; k)5(e', &,•), 
I 



(A4) 



the first of which is given in terms of a Gaunt-like coef- 
ficient, Q, and a projection coefficient onto a harmonic 
function, a^t- In the second equation, is multiplied 
by the LMTO structure constants, S, rescaled by the 
matching coefficient matrix, fi, that matches the atomic 
solutions inside the muffin tin to the Bessel and Neu- 
mann functions in the interstitial region (as well as their 
respective derivatives). The quantities A4t and Mht de- 
scribe the partial occupancies inside a muffin tin, and we 
may also use them to calculate, for example, the total 
orbital moment, 



M 



m-A4t(e, £, m, e, £' , m') 



e,e' ,t,m 



(A5) 



O t (e,e',£,£')8{£,£')5(m,m'), 



or, by summing away the energy set indices, we can get 
the total ^-charge, 



N e = M t (e,£,m,e',£' ,m') 

e,e',m 

O t (e,e',£,£')5(£,£')6{m,m'). 



(A6) 



Here we also need to account for O t (e, e', £, £'), the over- 
lap between the radial basis functions on type t. 



2. In the intersititial region 

The density in the interstitial is constructed from the 
k-dependent density matrix Pij(k) from Equation (|Aip 
and pseudobeLsis functions ^i(k), which are built with 
the constraints of matching the true basis functions in 
the interstitial and having a rapidly convergent Fourier 
representation. The ij) are, in practice, defined by their 
Fourier representation, and we have 

fi{r) = u(k, r) 
fc 

n(k,r) = Tr^V^(k,r)p. y (k)^](k,r) (A7) 

ij 

= Tr E F ^ ( k ' 9)}pij (k)^ 1 [$ (k, g)} 

where F~ l is the inverse Fourier transform. The pseudo- 
density n(r) is identical to the true density n(r) in the 
interstitial. The construction of ijj is described in Ref. 

n 

Because the density matrix is Hcrmitian, the density 
is accumulated for each k by 

fi(k, r) = Tr F^bPi (k, g)] ( £ 2 Pij (k)^ 1 [$ (k, g)\ 



j<i 



-pu^F-^l^gl 



(A8) 
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